Systematic Biology

Biologists often wish to explain why one trait is more common than its alternatives among species in a clade: is it due to asymmetries in diversification or character evolution processes? State-dependent speciation-extinction (SSE) models provide a flexible framework for jointly modeling such processes. Here, we derive mathematical results that define functional relationships (``rulesets") among rate parameters, establish a mapping between rulesets and evolutionary scenarios, and characterize how these scenarios give rise to different long-term tip state distributions. In particular, we show that multiple evolutionary scenarios can produce identical stationary distributions of tip state frequencies. Performing inference on simulated data from extant-only trees, we find that not all scenarios are equally favored for explaining datasets constructed solely from extant species. This bias could arise from the interplay between mixing rates, elapsed time, and the conditioning inherent in tree likelihoods. As a proof of concept, we apply our framework to the binary SSE (BiSSE) model across multiple Squamate clades and conduct a literature survey of previous BiSSE analyses on Angiosperm phylogeny. We show that, to some extent, the preference for certain scenarios observed in simulations is also reflected in analyses of empirical trees. This study provides a theoretical ground to directly link data-generating process and observed tip state frequencies under an SSE framework.

This study introduces a novel model that integrates the capacity of Potts models to capture structural dependencies with the computational efficiency afforded by Felsenstein's algorithm in phylogenetics. While the present work focuses on protein models, the PhyloPotts framework is equally applicable to alignments of discrete states, including RNA, DNA, and 3DI characters. Our findings indicate that, although phylogenetic estimates do not consistently exhibit reduced bias relative to standard models, the approach substantially enhances the accuracy of ancestral sequence reconstruction. The method is freely available via an open source implementation in the PhyloPotts package for BEAST 2 from https://github.com/rbouckaert/potts.

Species delimitation increasingly relies on multiple sources of evidence, but most frameworks are still limited to genetics and morphology or analyze data types separately and compare results qualitatively. We present delimSOM 2.0, an R package that leverages unsupervised machine learning through multilayer self-organizing maps. Any kind of input data, including genomic, phenotypic, ecological, and spatial, is analyzed jointly, retaining each data type as a separate layer and balancing their contributions through layer-specific weighting. Our package provides functions for preprocessing, replicate training, multivariate clustering, hyperparameter optimization, model-quality assessment, and variable- and layer-importance analyses. We evaluated our approach using two simulation sets and eight empirical studies spanning diverse taxa and data types. Overall, performance was strongest with larger sample sizes (>10 individuals per lineage) and modest numbers of lineages (K < 5), but unlike most species delimitation methods, it can reliably infer K = 1. Across two-population demographic simulations of genetic data, our method performed similarly to the well-established DAPC across gene-flow and sampling regimes, was more conservative than STRUCTURE and sNMF in inferring K = 2 and showed mixed support for K = 1 and K = 2 at intermediate differentiation levels and with gene flow. Crucially, our method is robust to missing data and variation in hyperparameters. Empirical examples showed that genetic data were typically most influential, but phenotypic, ecological, and spatial layers often contributed substantial complementary information. By combining diverse data types into a single unsupervised analysis, delimSOM 2.0 enables fully integrative species delimitation, reduces over-splitting, and advances efforts toward more objective and standardized taxonomy and speciation biology.

Early high disparity, where clades reach maximum morphological diversity early in their history, is often observed in the fossil record. Since there are currently no suitable methods to model these Early High Disparity (EHD) processes accurately, we introduce a new phylogenetic comparative method (PCM) to address this gap by combining two evolutionary modes--Ornstein-Uhlenbeck (OU) and Early Burst (EB)--leading to an initial surge in disparity until an early maximum, driven by the EB component, after which disparity declines due to OU constraints. We tested the EHD model on ichthyosaur skull length data and found it to be the best fit in nearly all cases compared to other modes. The median time of peak disparity (TEHD) is about 231.1 million years ago, during the Triassic, roughly 12.4% into the ichthyosaur lineage. To explore how environmental and geological factors might influence peak disparity, we performed changepoint analyses on sea level, sea surface temperature, rock area, and sedimentary sequence ages. While changepoints in temperature at 237.2 Ma and sequence boundaries at 238 Ma were closest to TEHD, no significant changepoints appeared to match TEHD.

Ubiquitous in phylogenomic and population-genomic introgression work, the D-statistic has fallen victim to its own success. Its reliability is compromised by two major confounding factors: among-lineage substitution rate heterogeneity, which generates spurious signals of introgression, and unsampled "ghost" lineages, which mislead the characterization of introgression by erroneously attributing ghost ancestry to extant species. Here, we developed D-plus, a maximum-likelihood tool that explicitly accounts for both confounding factors to unpack the sources of significant D-statistic signals. Simulations demonstrate that D-plus effectively mitigates false introgression calls under among-lineage rate variation. Meanwhile, it retains high sensitivity for genuine introgression signals and correctly assigns detected signals to introgression involving ghost lineages, or solely between extant non-sister lineages, or a combination of the two. Our D-plus re-analysis of phylogenomic data from Jaltomata (Solanaceae) overturns prior interpretations of significant D-statistic signals, which were conventionally regarded as evidence for recurrent non-sister introgression. Instead, our results favor alternative explanations invoking artefacts driven by lineage-specific rate heterogeneity as well as two ghost introgression events. Collectively, D-plus addresses major methodological limitations inherent to traditional D-statistic analyses and enables accurate detection and characterization of introgression events.

Motivation: Tree sets appear naturally when performing model based phylogenetic analyses. When tracking host-parasite co-evolution, testing sensitivity of different model assumptions, or even just checking convergence of individual runs of the same analysis, the question arises: are a given pair of tree sets the same? And if they are not the same, how do they differ? Results: We provide DensiTree 3, an interactive tool for inspecting pairs of sets of trees that allow visualising both tree sets using a single panel. Incorporation of CCDs provide better summary trees, taxon orderings and extra functionality like rogue detection, making it easier to interpret tree samples. Furthermore, a clade comparison graph allows inspection of differences in clade support and clade heights linked to the double-DensiTree panel. Together, these visualisations of tree sets make it easy to identify potential differences in the tree sets and interactively inspect these differences. Furthermore, recent advances in conditional clade distributions (CCDs) to represent tree distributions have been incorporated in DensiTree 3, and the most visible part of the user interface has been updated for a more modern look. Availability: DensiTree 3 is freely available from https://cs.auckland.ac.nz/~remco/DensiTree/. The program is licensed under GPL and source code is available via the link above.

Phylogenetic trees inferred from different genomic regions, or from organellar and nuclear genomes, frequently disagree owing to incomplete lineage sorting, hybridization, introgression or horizontal gene transfer. Such discordance is informative, yet the standard tool for comparing trees visually, the tanglegram, is restricted to two trees and is known to misrepresent their similarity. Here we introduce phylogenetic parallelograms, drawings in which any number of rooted trees are embedded into a common scaffold, a rooted phylogenetic network that displays all of them, and are drawn in parallel so that shared branches align while conflicting branches visibly diverge. On synthetic tree pairs, the complexity of a parallelogram tracks the topological distance between the trees, whereas the complexity of an optimized tanglegram does not. We apply the approach to genome-wide introgression in the Anopheles gambiae species complex and to organellar discordance in cats, grasses and Fagaceae. PhyloParallelograms is an open-source interactive application.

Speciation can be a slow process for long-lived organisms with large populations due to persistent incomplete lineage sorting, while hybridisation can further obscure species boundaries through admixture. In eucalypts, there is tremendous morphological and environmental variation with almost a thousand described species, subspecies, and hybrid taxa falling into multiple subgenera and intermediate sections. Multilocus molecular studies have generated concatenated species trees supporting deeper divergence between the Corymbia + Blakella + Angophora clade with the rest of eucalypts. However, the proportion of gene trees that recover the concatenated tree topology, including the taxonomic relationships among individual subgenera and sections, remains underappreciated. In this study, we assessed the predictive power of current eucalypt taxonomy - represented as a concatenated tree - in explaining topological variation across gene trees. To do this, we first extracted 1,187 BUSCO loci from short-read data of 701 samples representing roughly 500 described species using CAPTUS. We then inferred individual gene trees and a concatenated species tree using IQ-TREE2. Finally, we calculated the concordance factors for every described subgenus and section on the concatenated tree, and developed a taxonomic consistency test that tracks the taxonomic placement of each sample across gene trees. We show that, even though most taxonomic groups are monophyletic on the concatenated tree with 100% bootstrap support, several groups exhibit very low concordance, indicating that most gene trees do not recover their relationships as represented by the concatenated tree. These results highlight important considerations for current eucalypt taxonomic classification and provide a phylogenomic inference pipeline that can be applied not only to future eucalypt studies, but also to other evolutionarily complex groups of taxa.

Universal Single-Copy Orthologs (USCOs) are taxonomically-defined gene sets originally proposed to evaluate completeness in genome assemblies. Due to their single- or low-copy nature, they have also been proposed as ideal markers for phylogenomics. Here, we test this assumption using two vertebrate datasets, mammals and ray-finned fish, and their corresponding USCO reference marker sets, actinopterygii_odb10 and mammalia_odb10. Our results show that although most USCO markers are single-copy, 85% and 91% respectively, empirical analyses using a large number of taxa are able to identify a non-negligible proportion of paralogs in the actinopterygii_odb10 and mammalia_odb10 BUSCO gene sets. Gene phylogenies show that several of these paralogs derive from ancient whole-genome duplication events, particularly the teleost-specific whole-genome duplication, and also more recent duplications such as those in salmonids and cyprinids. Paralogs do not reflect the organismal phylogeny but rather the history of gene duplications, so retaining paralogs in phylogenomic datasets can lead to inaccurate phylogenomic hypotheses. To test the effect of ignoring paralogy in BUSCO gene sets, we analyze the mammal and ray-finned fish datasets using concatenated maximum likelihood and summary coalescent methods by randomly choosing one paralog per locus if these are present. The resulting phylogenomic trees are overall consistent with the reference phylogenies, although small discrepancies were found. Phylogeny-guided paralogy removal also changes the sequences retained for many loci, but whether these changes alter the inferred species trees remains to be tested.

The Asparagaceae is an economically and ecologically important angiosperm (monocot) clade, containing iconic garden plants such as asparagus, hyacinth and agave, and arid-adapted keystone species like Joshua tree and dragon blood tree. The approximately 2,500 species in the family are spread across seven subfamilies: Agavoideae, Aphyllanthoideae, Asparagoideae, Brodiaeoideae, Lomandroideae, Nolinoideae and Scilloideae. Relationships among the seven subfamilies are generally unresolved and the sources of gene tree-species tree discordance are largely unknown. Investigating the biology underlying such discordance requires a robust sampling of both taxa and molecular markers, due to the extremely short branches that form the backbone of the Asparagaceae phylogeny. Using targeted sequencing of 1,726 low-copy nuclear genes, and a nearly complete sampling of genera, we interrogate relationships among all seven Asparagaceae subfamilies and investigate the underlying sources of discordance across the phylogeny. Analyses using different phylogenomic approaches (multispecies coalescence vs. nuclear gene supermatrix vs. chloroplast gene supermatrix) each revealed extensive discordance across the Asparagaceae backbone, which may be associated with rapid lineage diversifications that occurred in Africa shortly following the origin of the family ~49 Ma. Specifically, we found that two branches, representing ancestors of Aphyllanthoideae and Nolinoideae, may be responsible for much of the backbone discordance. The placement of Aphyllanthoideae, represented by one extant taxon (Aphyllanthes monspeliensis), is controversial with discordant topologies across many studies. However, multispecies coalescence-based analyses from this study strongly support a sister relationship between Aphyllanthoideae and a larger Agavoideae + Scilloideae + Brodiaeoideae clade. According to gene tree quartet frequencies and ABBA-BABA tests, we cannot reject incomplete lineage sorting (ILS) as the main source of discordance for the placement of Aphyllanthoideae; our results leave open the possibility of an ancient reticulation involving an ancestor of Aphyllanthoideae and an extinct basal Asparagaceae lineage. Regarding Nolinoideae, we found that conflicting support for two possible bifurcating clades (Asparagoideae + Nolinoideae vs. Asparagoideae + Lomandroideae) was associated with skewed support for alternative gene tree topologies, suggesting that phenomena other than ILS may explain the observed discordance. A quartet-based ABBA-BABA analysis, used to test for signatures of introgression, revealed a putative ancient hybridization event between two divergent, ancestral lineages which may have contributed to the origin of the Nolinoideae clade in Africa ~45 Ma. In sum, this study elucidates the long-troublesome backbone of the Asparagaceae phylogeny and advances our understanding of the evolutionary and biogeographical origins of all extant subfamilies. More broadly, this study expands our understanding of the various sources of gene tree discordance that may arise from natural processes associated with rapid lineage radiations across the tree of life.

Chromosome-scale conservation of gene linkage (macrosynteny) provides phylogenetic characters largely independent of sequence data. Fusion of ancestral linkage groups, and especially the subsequent interspersion (mixing) of their genes, is argued to be effectively irreversible, and such synapomorphic events have been used to resolve deep, contentious nodes. I present Anvaya, a dependency-light, open-source toolkit of two Python programs with companion threshold-calibration tools which takes reciprocal-best-hit (RBH) synteny data directly to a phylogenetic analysis that codes fusion and mixing as discrete characters under explicit optimality criteria. The syntenybuilder module builds RBH orthology tables from genome assemblies, using either a designated reference genome or a multi-genome consensus procedure that reconstructs linkage groups de novo, requiring no curated reference set. The character-maps module codes every pairwise linkage-group co-occurrence as unfused, fused-clean or fused-and-mixed from gene order and a Fisher's exact test, then searches for the most parsimonious trees under an unordered (Fitch/Sankoff) and an ordered, irreversible (Camin-Sokal) cost model, optionally under a backbone constraint. Every tree is accompanied by ancestral states, synapomorphies mapped onto clades, consistency and retention indices, bootstrap support, and a permutation test of phylogenetic signal. Anvaya also detects whole-genome and segmental duplication from an intragenomic comparison and recovers the ohnologs that reciprocal best hit structurally drops, which can otherwise render a retained linkage group invisible. Outputs include vector figures and matrices in NEXUS and TNT/Hennig86 formats for verification in PAUP* and TNT. On synthetic data, consensus reconstruction resists a conflicting genome, rescue reverses a total false negative, and the diagnostics flag a deliberately homoplastic character. Across thirteen acarine and outgroup genomes, weakly supported associations concentrate in the most distant genome and in poorly assembled regions. In two duplicated chelicerates, ohnolog rescue recovered 12.5% (Limulus polyphemus) and 2.8% (Mastigoproctus giganteus) and restored three Limulus chromosomes that otherwise carried none.

Lineages that have undergone little apparent evolutionary change over long geological intervals are often described as "living fossils", with iconic examples including the coelacanth, tuatara, and Ginkgo. However, studies of morphological and molecular data have produced mixed evidence regarding whether living fossils have evolved more slowly than their relatives. We tested for reduced evolutionary rates using a sister-pairs approach, analysing a comprehensive set of morphological and molecular data for 45 living-fossil lineages and their close relatives across the Tree of Life. Living-fossil lineages showed substantially reduced rates of morphological evolution, whereas evidence for reduced molecular rates was weaker. To test whether this pattern could be explained by species poverty, we also analysed other depauperate lineages (species-poor clades contrasted with species-rich sister clades) and found a weak tendency towards reduced molecular rates but without a corresponding slowdown in morphological evolution. Our study reveals that the living-fossil concept captures a quantifiable pattern of prolonged morphological conservatism, but also suggests that living fossils lack other shared biological or life-history characteristics that might unite them as a group.

Ancestral recombination graphs (ARGs) represent the coalescent history at every locus in the genome, subject to mutation and recombination. Modelling viruses using the coalescent with recombination remains less explored, due to the potential for complex, multi-scale dynamics. Using perturbative methods, we establish mathematical bounds that describe the balance between inter-host and within-host coalescence, and demonstrate that the highly-recombinant Epstein-Barr Virus (EBV) coalesces predominantly at the inter-host level. We then modify SINGER, a tool from diploid genetics, for haploid ARG inference with fine-grained variable recombination rates. We provide insight into the paradox of high recombination rates and slow linkage disequilibrium decay, as ARG-based methods potentially allow us to disentangle epistatic effects from shared ancestry, and illustrate this using the coalescent topologies of EBV's host cell-entry glycoproteins. Our findings demonstrate that ARG-based methods recover gene-level selective structure, and provide a general method for analysing virus phylodynamics in the limit of small quantities.

The controversy over Elateriformia relationships hampers understanding of the evolution of disparate soft-bodied and fully sclerotized phenotypes, bioluminescence, and the ontogenetic reprogramming that leads to larviform females. Elateriformia contains >50,000 species in 34 families, including fireflies, jewel, and click beetles. We sample 31 families and 195 ingroup taxa, each with up to 4,224 orthologs. Through explicit hypothesis testing, coalescent-based approaches, and evaluation of dataset properties, we resolve backbone relationships. We reconfirm the monophyly of Byrrhoidea sensu lato, rejecting the recent resurrection of Dryopoidea, and verify the reciprocally monophyletic lampyroid and elaterid clades. We characterize short internodes and gene-tree discordance, both of which are indicative of ancient rapid radiations. Counterintuitively, signals from moderately variable genes are crucial for resolving certain splits, suggesting that even subtle rate differences distinguish phylogenetically informative loci from overly conserved, less informative markers. Time-calibrated phylogenies indicate the early origin of Elateriformia (281-287 Mya), with superfamilies diverging at the Permian-Triassic boundary and diversifying throughout the Jurassic. We identify 15 independent origins of soft-bodiedness and neotenic modifications spanning the Mid Jurassic to the present. Our hypothesis provides a foundation for investigating the genetic mechanisms underlying ontogenetic reprogramming and bioluminescence.

Trade-offs are commonly assumed when the statistical parameters used to define biological traits are negatively correlated. However, parameter correlations can emerge from the geometry of the model, rather than underlying biological signal. We characterize this phenomenon generally, and develop a framework to parse biological signal from statistical artifact in parameter correlations. Using the widely reported trade-off between thermal sensitivity and static upper critical thermal limit as a case study, we show that an extreme correlation between these parameters is expected even in the absence of a real trade-off. However, our method reveals that there is, in fact, a biological correlation underneath this spurious artifact. We confirm this finding, as independent genome-wide association analysis reveals genetic variation among trade-off phenotypes, but none corresponding to net thermal tolerance, suggesting that adaptive capacity to a warming climate may be severely limited. This particular trade-off is therefore biologically supported, illustrating how biological signal can be obscured by estimation artifact, thus evading correlation-based detection. We broaden our analysis, showing that artifact is present across biological systems and model types, although its manifestation may vary. Analytical results and simulations reveal how sampling design, model choice and parameterization govern the artifact, providing practical guidance for experimental planning. Our findings identify both a problem and a solution: correlations among model parameters cannot be interpreted biologically without first accounting for the structure imposed by the model and data, but genuine biological relationships can still be recovered when that structure is explicitly modeled.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Three-dimensional imaging is rapidly increasing the scale of biological morphometric datasets, but establishing automated anatomical correspondence remains a major bottleneck. Conventional point-cloud registration methods transfer landmarks efficiently, yet poor initialization or high morphological variation can produce geometrically plausible but biologically implausible deformations. Here, we developed a statistical-shape-model (SSM)-regularized coherent point drift (CPD) estimator for automated morphometric landmark transfer. It replaces the unconstrained non-rigid CPD update with optimization in SSM coefficient space and regularizes coefficients by inverse population covariance, so deformations along low-variance biological modes are penalized more strongly. We compared pose-matched conventional CPD, SSM-CPD, and SSM-CPD with refinement in leave-one-out atlas experiments on laboratory mouse skulls (N=249), primate molars (N=115), and squamate vertebrae (N=254) under aligned, mildly perturbed, and severely perturbed initialization. We evaluated these formulations in terms of landmark accuracy, failure-tail behavior, runtime, and preservation of morphometric structure. Under severe perturbation, we also evaluated an integrated recovery procedure that optimizes the starting template, recovers global pose, and then applies SSM-CPD with refinement. Severe misalignment produced widespread failures in pose-matched conventional CPD, whereas SSM regularization reduced median error from 13 - 21% to approximately 1 - 3% of the target bounding-box diagonal. This improvement was highly consistent for skulls and molars. In the morphologically complex vertebral dataset, however, low median error concealed a substantial minority of catastrophic failures (28% of vertebrae above a 6% error threshold), which integrated recovery reduced to 4%. Across datasets, automated landmarks preserved the major axes and pairwise relationships of morphometric shape space, and local CPD refinement improved covariance preservation. Embedding learned biological variation directly into registration broadens the range of starting conditions and anatomies for which automated landmark transfer remains reliable. Our results also show that failure tails and downstream morphometric fidelity must be assessed alongside median accuracy when automated landmarking methods are evaluated.

Flowering plants vary in sex allocation (resource division between male and female function) and in selfing rate. These traits should evolve together: selfing changes fitness returns from pollen and ovules, whereas sex allocation changes the pollen available for selfing and export. Yet most models allow only one trait to evolve, leaving their evolutionary feedback unresolved. Here we model the joint evolution of sex allocation, pollen export and maternal control of competition between self- and cross-pollen. We show that selfing and sex allocation strongly influence one another, as selection favours associations between sex allocation and mating-system traits. When maternal control is fixed to recover classical mass-action, evolution can lead to the stable coexistence of selfing hermaphrodites and males specialised in pollen export. This androdioecy persists because hermaphrodites cannot prevent fertilisation by males. When maternal control of pollen competition also evolves, we find that where a mixed-mating equilibrium exists, selection again favours the emergence of androdioecy. As hermaphrodites evolve to favour their own pollen, however, they progressively exclude pollen from males, driving them extinct. The population may then return to mixed mating and repeat the process, or evolve complete selfing. We discuss these results with reference to known variation in the mating and sexual systems of plants.

Metabolic rate governs the pace of life through the allocation and expenditure of energy, yet metabolic rate data are unavailable for most species. External morphological traits, such as gill slits, have been used to infer respiratory physiology and metabolism in sharks, yet, whether that pattern is generalizable across rays is unknown. In rays (Batoidea), spiracle length seems adaptive to oxygen uptake capacity, environmental oxygen availability, and activity level. Thus, we evaluate the allometry of external respiratory morphology (spiracles and gill slits) in rays using 550 museum specimens representing 88 species. We test whether spiracle length and gill slit height grow isometrically with body size (total length or disc width) within and across species. Overall, we find spiracles and gill slits are morphologically diverse, varying widely in shape within and between families. Spiracle length and gill slit height scale isometrically and are a consistent proportion of body size, despite such high morphological variability. While that proportion is consistent through ontogeny, each species proportion differs (e.g. gill slit height ranges from 0.93% of total length in Spotted Legskate to 2.93% in Pacific Electric Ray). Of the three swimming styles, pectoral oscillating species (Mobulid and Eagle Rays) may have steeper scaling of spiracle length with gill slit height - consistent with greater respiratory capacity in species with more active lifestyles. Our findings are consistent with literature hypothesizing external respiratory morphology (spiracles and gill slits) is tightly related to metabolic rate and activity level, and demonstrate the value of natural history collections in ecological research. Importantly, the isometry of both respiratory traits means that a single measure of any individual, irrespective of size, will be representative of the species thus providing a basis for expanding our understanding of metabolic morphology, enabling future work on life history, extinction risk, and ecological function.

Dispersal and hybridization can reinforce one another in a feedback loop when species that disperse farther are more likely to encounter and mate with close relatives, resulting in gene flow that can reintroduce dispersal-enhancing alleles in their descendants. This dynamic could help explain why complex dispersal-related traits, such as flight in insects, appear to evolve repeatedly within a clade, but this idea has rarely been tested with genomic data. North American Gryllus field crickets, in which flight capability appears to have been regained from flightless ancestors at least nine independent times, offer a powerful system to do so. Using a genome-wide dataset, we test three predictions about the origins of phylogenetic discordance in Gryllus and the role of dispersal-related traits. First, we find that both incomplete lineage sorting (ILS) and hybridization contribute to extensive gene tree conflict observed across the Gryllus phylogeny, but several conflicting patterns cannot be explained by ILS alone. Multiple lines of evidence reveal introgression at both deep and recent nodes, producing a mosaic evolutionary history with reticulation across many lineages. Second, hybridization is disproportionately concentrated among flight-capable lineages, supporting the idea that dispersal ability increases opportunities for interspecific gene flow. Third, introgression from flight-capable lineages could help explain repeated trait reversal across the phylogeny, suggesting that allele flow may have contributed to the re-emergence of flight in descendants of flightless ancestors. These results are consistent with dispersal-hybridization feedback shaping both the reticulated phylogenetic history of Gryllus and the repeated emergence of flight. More broadly, this species group reveals how models of character evolution built on bifurcating trees may overestimate the number of independent trait origins when reticulate processes are ignored.